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ABSTRACT 

We study the fate of internal gravity waves approaching the centre of an initially non- 
rotating solar-type star, primarily using two-dimensional numerical simulations based 
on a cylindrical model. A train of internal gravity waves is excited by tidal forcing 
at the interface between the convection and radiation zones of such a star. We derive 
a Boussinesq-type model of the central region of a star and find a nonlinear wave 
solution that is steady in the frame rotating with the angular pattern speed of the 
tidal forcing. We then use spectral methods to integrate the equations numerically, 
with the aim of studying at what amplitude the wave is subject to instabilities. These 
instabilities are found to lead to wave breaking whenever the amplitude exceeds a 
critical value. Below this critical value, the wave reflects perfectly from the centre of 
the star. Wave breaking leads to mean flow acceleration, which corresponds to a spin 
up of the central region of the star, and the formation of a critical layer, which acts 
as an absorbing barrier for subsequent ingoing waves. As these waves continue to be 
absorbed near the critical layer, the star is spun up from the inside out. 

Our results point to an important amplitude dependence of the (modified) tidal 
quality factor Q', since nonlinear effects are responsible for dissipation at the centre 
of the star. If the amplitude of the tidal forcing exceeds the critical amplitude for 

wave breaking to occur, then this mechanism produces efficient dissipation over a 

1 

continuous range of tidal frequencies. This requires ^ jjy^ ( i Tay ) ^ ^ ^'^"^ ^ planet 
of mass TOp in an orbit of period P around the current Sun, neglecting stellar rotation. 
However, this criterion depends strongly on the strength of the stable stratification 
at the centre of the star, and so it depends on stellar mass and main-sequence age. 

If breaking occurs, we find Q' w 1.5 x 10^ (r£ij;) ^ ' current Sun. This varies 

by no more than a factor of 5 throughout the range of solar-type stars with masses 
between 0.5 — I.IM0, for fixed orbital parameters. This estimate of Q' is therefore 
quite robust, and can be reasonably considered to apply to all solar- type main-sequence 
stars, if this mechanism operates. The strong frequency dependence of the resulting 
dissipation means that this effect could be very important in determining the fate 
of close-in giant planets around G and K stars. We predict fewer giant planets with 
orbital periods of less than about 2 days around such stars if they cause breaking at 
the centre, due to the efficiency of this process. 

Even if the waves are of too low amplitude to initiate breaking, radiative damping 
could, in principle, lead to a gradual spin-up of the stellar centre and to the formation 
of a critical layer. This process could provide efficient tidal dissipation in solar-type 
stars perturbed by less massive companions, but it may be prevented by effects that 
resist the development of differential rotation. 

These mechanisms would, however, be ineffective in stars with a convective core, 
such as WASP-18, WASP-12 and OGLE-TR-56, perhaps partly explaining the survival 
of their close planetary companions. 

Key words: planetary systems - stars: rotation - binaries: close - hydrodynamics - 
waves - instabilities 
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1 INTRODUCTION 

Tidal interactions are thought to be important in determin- 
ing the fate of short-period extrasolar planets and the spins 
of their host stars, as well as in the sychronization and cir- 
cularization of close binary stars. The extent of spin-orbit 
evolution that results from tides depends on the dissipa- 
tive properties of the bodies involved in the interaction. It 
is standard to parametrize our uncertainties in the mecha- 
nisms of tidal dissipation in each body by defining a tidal 
quality factor Q, which is an inverse measure of the dissipa- 
tion. This is usually defined to be proportional to the ratio 
of the maximum energy stored in a tidal oscillation (Eq) to 
the energy dissipated over one cycle, i.e. 



Q = 2nEo 



-Edt 



(1) 



and we find it convenient to use the variant Q' = where 
k is the second-order potential Love number of the body. 
This combination always appears together in the evolution- 
ary equation^ 

Q' is difficult to calculate from first principles in fluid 
bodies, and uncertainties in the mechanisms of tidal dissi- 
pation remain. The tidal disturbance can generally be de- 
composed into two parts: an equilibrium tide and a dynam- 
ical tide. The equilibrium tide is the quasi-hydrostatic el- 
lipsoidal tidal bulge. In the frame corotating with the fluid, 
the time-dependence of the equilibrium tide is dissipated 
through its interaction with turbulent convection, though 
the damping rate is uncertain, particularly when the con- 
vective time exceeds the tidal period ( Zahn|1966 Goldreich 



|fc Nicholson"lT977l [Goodman fc Oh|1997[|Penev et al.|2007| ) 

The dynamical tide consists of internal waves that are ex- 
cited by low-frequency tidal forcing, and has received much 
recent interest with regard to its possible contribution to 
Q' ("Witte & Savonijc 20021 lOgilvie & Lin 



2004 



hereafter 



OL04; ,Wu,.2005, Papaloizou fc Ivanov||2005| [ Ivanov fc P~- 
paloizou||2007| lOgilvie fc Lin||2007[ Hereafter OL07; |Goo~ 



man fc Lackner||2009 1. This is because if these waves have 
short wavelength, then they are more easily damped than 
the large-scale equilibrium tide by radi ative diffusion (|Zahn| 
1975[|Zahn|1977[ ), con vective viscosity ([Terquem et al.|1998[ ) 
or nonlinear breaking ( Goodman fc Dickson|19^ hereafter 
GD98). 

The tidal frequency is typically much lower than the dy- 
namical frequency of the body, so the relevant internal waves 
must be approximately incompressible, restored not by pres- 
sure but by buoyancy or rotation. OL04 found that the dissi- 
pation of tidally excited inertial waves, whose restoring force 
is the Coriolis force, can contribute signiflcantly to the dis- 
sipation rate in a giant planet, whose interior is mostly con- 
vective (see also Wu||2005 and Ivanov fc Papaloizou||2007 |. 
These waves can also contribute to the dissipation rate in 
convective regions of stars (OL07). They are excited by tidal 
forcing of frequency to, if this is less than the Coriolis fre- 
quency (2n), and this is true for many astrophysically rele- 
vant circumstances. However, these waves are not excited if 
the tidal frequency exceeds the Coriolis frequency, so this 
process is then not effective at dissipating the tide, and 



contributing to Q' . In this work we concentrate on waves 
that have buoyancy as the restoring force, and which are 
commonly referred to as internal gravity waves (hereafter 
IGWs^Q 

IGWs have been proposed to account for the efficient 
tidal dissipation that has been inferred from the circularisa- 
tion of early- type binary stars ( Zahn|1975[ Zahn|1977| Zahn 
'2008; Savonijc fc Papaloizou 1983; Papaloizou & Savonije 



1985, ,Savonije et al. 1995; Savonije fc Papaloizou|1997[ Pa- 



paloizou fc Savonije 19971, which are massive enough to 



have a convective core and an exterior radiative envelope. 
In these stars, IGWs are excited near the boundary be- 
tween these two regions, where the buoyancy frequency (or 
Brunt- Vaisala frequency, see § [2]| matches the tidal forcing 
frequency. These propagate outwards into the stably strati- 
fied radiation zone, towards the surface, where they are fully 
or partially damped by radiative diffusion. In this picture, 
these stars are tidally synchronized from the outside in, since 
angular momentum is deposited in the regions of the star 
where these waves damp ( Goldreich fc Nicholson|1989 a.,b). 

The above model only works for stars with an exterior 
radiation zone, which is unlike that of the Sun and other 
stars of solar type, which have radiative cores and convective 
envelopes. It is of particular interest to study the efficiency of 
tidal dissipation in these stars, since many have been found 
to harbour close-in planets, whose survival is determined by 
the stellar Q' . This is because a planet with an orbital pe- 
riod shorter than the stellar spin period is subject to tidally 
induced orbital decay, with an inspiral time that depends on 



dissipation in the star (e.g. Barker fc Ogilvie 2009 OL07) 
In these stars, a train of IGWs are again excited at the in- 
terface between the convective and radiative zones, but here 
they propagate towards the stellar centre. If they can coher- 
ently reflect from the centre, global standing modes can form 
in the radiation zone. In this case, tidal dissipation is effi- 
cient only when the tidal frequency matches that of a global 
standing mode (which are commonly referred to as g-modes) 
( |Terquem et"ar|p98l [Savonije fc Witte||2002[ ). This would 
not contribute appreciably to Q because the system would 
evolve rapidly through these resonances, unless resonance 
locking occurs ( [Witte fc Savonije|p99| [Witte fc bavomje| 
20011. On the other hand, if these waves do not reflect co- 



herently from the centre, and are either strongly dissipated 
there, or are reflected with a perturbed phase, then efficient 
dissipation is possible over a broad range of tidal frequencies 
(GD98; OL07). The extent of nonlinearity in the waves near 
the centre is likely to be the factor that determines whether 
these waves reflect coherently, and this is controlled by the 
amplitude and frequency of the tidal forcing, as well as the 
properties of the stellar centre (OL07). 

In this paper we study the problem of IGWs ap- 
proaching the centre of a solar-type star, primarily using 
two-dimensional numerical simulations. We first derive a 
Boussinesq-type system of equations appropriate for the 
stellar centre, which are ideal for integrating numerically 
using spectral methods. An exact solution for tidally forced 
waves is derived, and some of its properties are discussed. 
Our numerical set-up is described and results are presented 



^ Q' reduces to Q for a homogeneous fluid body, where fc = | . 



^ Though note that they are often referred to as g-modes, or 
simply gravity waves, in the literature. 
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for both linear and nonlinear forcing amplitudes, including 
an analysis of the reflection coefficient and a study of the 
growth of different azimuthal wavenumbers in the distur- 
bance. This is followed by a discussion of the results, espe- 
cially their relevance to Q' for solar-type stars, and to the 
survival of close-in giant planets in orbit around such stars. 



2 INTERNAL GRAVITY WAVES: 

ELEMENTARY PROPERTIES, WAVE 
BREAKING AND CRITICAL LAYERS 

IGWs are a family of dispersive waves that are ubiquitous 
in nature. They propagate in any fluid with a stable density 
stratification, due to the restoring force of buoyancy. Their 
influence can be observed in the oceans and atmosphere of 
the Earth on a range of spatial and temporal scales, from 
the visual undulations of striated cloud structures, to the 
complex interplay between these waves and shearing flows, 
which produces the large-scale Quasi-Biennial Oscillation 
in the equatorial stratosphere. It is widely recognised that 
IGWs play a prominent role in the transport of energy and 
angular momentum in geophysical and astrophysical flows 
Mclntyro^ 2000; Rogers & Glatzmaierl[2006 



(Biihler 2009 



ILJ . — :^ II 1 I II 
19991. IGWs are also thought to be important 



in stably stratified radiation zones of stars. When excited 
by turbulent convection, they were at one stage put forward 
as potential explanations for maintaining the solid body ro- 
tation of the radiative interior of the Sun ( |Schatzman|1993 



Zahn et al. 19971. However, it was pointed out that the 



"antidiffusive" nature of IGWs tends to enhance local shear 
rather than reduce it (Gough & Mclntyro 1998). IGWs are 



still thought to produce angular velocity variations in the 
radiation zone ( Rogers fc Glatzmaier|2d06 1. They have also 
been invoked to explain the Li depletion problem in F-stars 
( [Garcia Lopez &: Spruit|1991[ ), affecting solar neutrino pro- 
duction ( Press||1981 1, and possibly having an effect on the 
solar cycle ( Kumar et al.|19"99 L 

Observations of oscillations on the solar surface are able 
to provide information about the interior properties of the 



Sun ( Christensen-Dalsgaard 20021. IGWs in the radiative 
interior of the Sun can form global standing modes, com- 
monly referred to as (j-modes, if their frequency matches 
that of a free mode of oscillation. These are known to have 
their amplitude largest close to the centre (see solution in 
§ [sj and would therefore seem ideal probes of the deep inte- 
rior. Unfortunately for observers, the standing jr-modes are 
effectively trapped in the radiative interior, where the strat- 
ification is stable, and are evanescent in the convection zone, 
and so are unlikely to be visible at the solar surface. Never- 
theless, modes of sufflciently low degree, with high enough 
amplitude, may have already been observed at the surface 



by Garcia et al. ( 2007 1, though it must be noted that thus far 



there is no undisputed evidence for observations of g-modes 
( Appourchaux et al.|2009 l. 

The frequencies of the largely incompressible internal 
waves lie in ranges controlled by the buoyancy frequency 
(or Brunt- Vaisala frequency) A'' and the Coriolis frequency 
2Q,. The square of the buoyancy frequency in a spherically 
symmetric star is defined by 



N^{r) = 



1 dp 
p dr 



1 d\np d\np 
Fi dr dr 



(2) 



where p,p are the density and pressure, and Fi = ^§j|-^j 

(at constant specific entropy s). The local dispersion rela- 
tion for linear noncompressive internal waves in a fiuid body 
rotating with angular velocity is 



up — N'^ sir? a + AT? cos^ /3, 



(3) 



where a is the angle between the wavevector k and the grav- 
itational acceleration g, and jS is the angle between k and 17. 
The frequency of these waves is independent of wavelength 
(in the absence of viscosity or thermal conduction) , and only 
depends on the direction of the wavevector. This is differ- 
ent from waves whose restoring force is due to compress- 
ibility, whose frequency is inversely proportional to wave- 
length. When A'^ = 0, Eq|4] describes inertial waves, which 
have frequencies in the range (0,2f2). If the body is non- 
rotating (11 = 0), then these waves are IGWs, and possess 
frequencies in the range (0, N). In the presence of nonzero 
and A^, these waves are intermediate between inertial waves 
and IGWs and are referred to as inertia-gravity waves. For 
waves in a spherical star, at a given latitude there is a min- 
imum frequency for inertia-gravity wave propagation. Near 
the equator, waves can propagate with arbitrarily low fre- 
quency. From here on we neglect the bulk rotation, and as- 
sume that f2 = 0, i.e., we consider only IGWs. The local 
dispersion relation for IGWs can be rewritten 



(4) 



where kh and kr are the horizontal and radial wavenumbers. 

The phase and group velocities of these waves can be 
calculated from Eq. |4] to give 

_ iu k _ Nkh 
'''' ^ |k| |k| " (fc2_^fc2)| 

- ■ ^^'^ {krGr + kneh) 

Nkr 



{kre,. + kneh) 



k'^ 
_N 



2 {^h^r kr^h) 5 



{kh^r krGh 



(5) 

(6) 
(7) 



(8) 



in the tidally relevant limit that the radial wavelength of 
the waves is much shorter than the horizontal wavelength, 
i.e., kh <C kr (which is true except near turning points or 
within the last few wavelengths from the centre of a star). 



In this limit, €„ ■ = —N-fy 



-Cr, ■ Br, i.e., the radial wave 



pattern moves in the opposite direction to the radial energy 
fiux. Since u) is independent of |k|, Cg • k = 0, meaning that 
the energy in IGWs propagates along surfaces of constant 
phase. 

For waves on a non-zero horizontal background shear 
fiow U, Eq. |4] still applies if the Richardson number Ri — 
|9u/9rp ^' if ^6 replace uj by the Doppler-shifted fre- 
quency u), and similarly for the phase and group velocity of 
the waves 



kU, 



U, 



u. 



(9) 



We now have the possibility of the wave frequency being 
Doppler-shifted upwards to A'^, in which case Cg ■ reverses, 
giving total internal refiection ( Mclntyre||2000 1 . 

The other extreme, that d) — >■ 0, occurs when the hor- 
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izontal velocity in the shear matches that of the horizon- 
tal phase velocity. This occurs at a so-called critical layer, 
which is defined as the layer at which the wavelength of the 
waves would be Doppler-shifted to zero, if they were ever 
to reach it. Note though, that Cg ■ as a) — >■ 0, so 

in linear theory the waves never reach the critical layer in 
a finite time. Early work on IGWs in a background shear, 
including a study of critical layers, can be found in IBooker] 
|fc Bretherton] ( |1967| ) and ,Hazel| ( |1967| ). They find that an 
IGW propagating through a critical layer is attenuated by 
a factor ~ exp (^-2-K{Ri - 1/4)^/^^ . If iii > 1/4, the wave 
is fully absorbed, and irreversibly transfers its energy to the 
mean flow. However, it must be noted that at the critical 
layer, linear theory predicts that the wave steepness 



max(u ■ e^j 

Cp ' 



(10) 



where u is the velocity perturbation, of these waves goes to 
infinity, as the Doppler-shifted horizontal phase speed Cp • 
goes to zero, i.e., the waves become strongly nonlinear at 



the critical layer (which Mclntyre 2000 refers to as linear 



theory predicting its own breakdown). Since wave breaking 
is expected to occur whenever s > 1, the waves are likely 
to break before they reach the critical layer. Nonlinear ef- 
fects have been studied in early simulations by Winters &' 



D'Asaro ( 19941, who find that the of initial wave energy on 



encountering a critical layer, roughly one third reflects, one 
third results in mean flow acceleration and the remainder 
cascades to small scales where it is dissipated. This implies 
that wave absorption by the mean flow need not be com- 
plete, in contrast with the prediction from linear theory. 
This is discussed further in § |9.2| in relation to the results of 
our simulations. 

Wave breaking is deflned as a wave-induced process that 
leads to the rapid and irreversible deformation of otherwise 
wavy material contours ( Mclntyre|2000 ), and it leads to the 
production of turbulence and irreversible energy dissipation. 
The breaking process results from the growth of an insta- 
bility upon a basic state composed of a wave with s > 1. 
The susceptibility of a wave to breaking can be enhanced 
by resonant triad interactions, in which a primary wave res- 
onantly interacts with a pair of low-amplitude secondary 
waves. This process transfers energy to the secondary waves, 
whose steepness can then grow beyond the critical value re- 
quired for breaking to occur, even though the primary steep- 



ness may not be sufficient for breaking on its own (Staquet 
fc Sommeria||20d2 I. Previous work has shown that the pro- 



cess leading to the wave steepening is two-dimensional, but 



that breaking is a three-dimensional process ( Klostermeyer 
1991|Winters fc D'Asaro 19941. Nevertheless, the mecha- 



nisms responsible for breaking and the final outcome of the 
breaking process are likely to be similar in 2D. 

Here we are interested in studying what happens when 
IGWs excited by tidal forcing approach the centre of a star 
with an inner radiation zone, i.e. G-type stars such as the 
Sun, which do not possess a convective core. We are primar- 
ily concerned with the tidal dissipation efficiency for solar- 
type stars that may result from breaking of these waves near 
the centre. Throughout the rest of this paper we restrict our 
problem to 2D, and postpone study of any 3D effects. We 
now describe our basic problem in more detail. 



3 BASIC DESCRIPTION OF THE PROBLEM 
3.1 Tidal potential 

Consider a star and a planet in a mutual Keplerian orbit 
(though we make no assumptions about the relative masses 
at this stage). The tidal potential experienced by the star 
can be written as a sum of rigidly rotating spherical har- 
monics (e.g. OL04). For the simplest case of a planet on a 
circular orbit, that is coplanar with the equatorial plane of 
the star, we can consider a two-dimensional restriction of 
the problem to this plane. This allows us to write the time- 
dependent part of the quadrupolar {I = 2) tidal potential in 
the equatorial plane of the star as 



*(r, <;/),*) = 



3 Gm 



r cos{2(f> ~ cut), 



(11) 



in spherical polar coordinates (r, 9, cj>) with origin at the cen- 
tre of the star, in the plane 6 — 7r/2. Here rrip is the planet 
mass, a is the orbital semi-major axis, and n > is the 
mean motion. The relevant tidal frequency is = 2n — 2fi 
for a star rotating with angular velocity $1. From here on 
we assume that the star is non-rotating, i.e. we assume that 
Q = 0, which is a reasonable assumption if the star is spin- 
ning much slower than the orbit. This is justified for the 
problem of a close-in gas giant planet with a several-day 
orbital period around a solar-type star that has been spun 
down by magnetic braking for the duration that it has spent 
on the main sequence, to rotate with a spin period of several 
tens of days. 

We consider a restriction of the full three-dimensional 
problem, in which the tidal potential is composed of many 
different spherical harmonic components for an orbit of arbi- 
trary eccentricity and inclination, to instead consider a sim- 
plified two dimensional model of a star, forced by this single 
component of the tide. The main motivation for this deci- 
sion is simplicity, since this is a first attack on the problem, 
though it is likely that this is a reasonable approximation. 
Extensions can be studied in the future. 



3.2 Central regions of a star 

The buoyancy frequency is real and comparable with the 
dynamical frequency of the star 



_ [ Gm, 



(12) 



throughout the bulk of the radiation zone, where m^, and 
Ri, are the stellar mass and radius, respectively. In Fig. [l] we 
plot A'^ normalised to cddyn in the radiation zone for Model S 
of the current Sun ( |Christensen-Dalsgaard et al. 11996] ). For 
our problem, the tidal frequencies of interest lj ^ N, which 
implies kr ^ 27r/r, except near the centre. IGWs are excited 
at the top of the radiation zone by a combination of tidal 
forcing in that region, together with the pressure of inertial 
waves acting at the interface if < 2Q (OL07). It is within 
this transition region that N increases linearly with distance 
into the radiation zone, so there is a point at which A ~ li, 
and IGWs are efficiently excited. These propagate towards 
the centre with radial wavelengths Ar < 10~^ — 10~^Rq, for 
typical tidal frequencies. 

Expanding the standard equations of stellar structure 
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Figure 1. Top: Buoyancy frequency N normalised to dynamical 
frequency ujiiy„ versus radius, based on Model S of the current 
Sun. Also plotted is a frequency 27r/lday ~ 0.1cj;jy„, correspond- 
ing with the orbital frequency for a one day orbit - note that 
this is half the tidal frequency, so aj <C Af throughout the bulk 
of the radiation zone. Only near the centre and at the interface 
between the convection and radiation zones, does li ~ Af. Bot- 
tom: Coefficient C in the expansion N Cr near the centre of 
the Sun normalised to its current value versus main-sequence age 
for a sequence of solar models that pass through Model S of the 
current Sun. The stratification steepens as the star evolves. 



about r = 0, we obtain the density stratification 

p(r) = po + P2r-'+0(/). (13) 

For sufficiently small r, <; is linear in r, and A'^ = Cr is 
also linear in r, where C is a constant that varies with 
stellar model and main-sequence age. For the current Sun, 
C = Cq ^ 8.0 X 10""m"^s~^ This is valid throughout 
only the inner < 3% of the Sun. However, even with the 
largest expected radial wavelength, this region contains mul- 
tiple wavelengths. In Fig. [l] we plot the variation in C ver- 
sus main-sequence age, normalised to its value in the cur- 
rent Sun, from a sequence of solar models that pass through 
Model S. The increase in C with age is due both to the in- 
creasing central condensation, and the build-up of a gradient 
in the hydrogen abundance (there is a small drop around 8 
Gyr when hydrogen is nearly used up and the contribution 
from the composition gradient decreases, after which the 
central density increases rapidly). 



4 DERIVATION OF A BOUSSINESQ-TYPE 
SYSTEM OF EQUATIONS 

The ideal compressible fluid equations in 2D plane polar 
(r, (j)) coordinates are 



DUr — ~-drP — 9r<&, 

r p 



UrU^ 



pr 



-d4,p 9,/,$, 



Dp + P 



-drirUr) + -dd,us 
r r 



0, 



Dp - '^Dp - 



0, 



D — dt + Urdr + —( 

r 

V^-l- = 47rGp. 



(14) 
(15) 

(16) 

(17) 

(18) 
(19) 



The basic state is static and circularly symmetric. Near the 
centre, we pose the expansion 



p = po + p2r^ + PAT*' + 0{r^), 



(20) 



and similarly for p and The coefficients are related by 
the condition for the background to be in hydrostatic equi- 
librium, together with Poisson's equation, giving 



P2 = -po4'2, Pi = -po$4 + 



P2p2 



$2 = TfGpo, '1>4 



nGp2 



(21) 
(22) 



and so on for the coefficients of higher-order terms. 

We are interested in a region where r/i?* = 0(e), where 
e ^ 1, so let r — ex. Now introduce a slow time r = et, then 
the solution including the basic state and a slow nonlinear 
density perturbation has the form 

p = po + e^P2X^ + e^pAx"^ + ... + <? P2{^, 0, t) + ... (23) 

V ' ^ ^ ' 

basic state nonlinear perturbation 

p = po + ep^x' + ep^x^ + ... + e^p4(a:, 0, r) + ... ,(24) 



basic state 

basic state 
Ur = e'^Ur2{x,(j>,T) + 



nonlinear perturbation 

= e^u^2{x, (j>, t) + 



nonlinear perturbation 

.+ e%'^{x,^,r) + .{2^) 

' ^ V ' 

nonlinear perturbation 

(26) 



(27) 



nonlinear perturbation 

In the Boussinesq approximation, the fractional pressure 
perturbation is small compared with the fractional density 
perturbation because this is a low-frequency perturbation, 
with u) <^ ijJdyn, for which acoustic effects are negligible. For 
such short-wavelength perturbations, the gravitational per- 
turbation is also small ( Christensen-Dalsgaard|2002 1 . These 
two points explain the absence of terms in the above solu- 
tions proportional to for p and $ in the nonlinear pertur- 
bation. This corresponds to looking for small density con- 
strasts and weak material accelerations compared with grav- 
ity. 

Substituting these expansions into the basic equations, 
and subtracting terms that arise only in the basic state, we 
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obtain at leading order 

DlUr2 



DlU^2 + 

1 



X 

Ur2U4,2 



1 , 2xp2p'2 o ^/ 
— OxVi H 5 Ox'^i, 

Po 



Pi 



Pox 



-d<t>p'i 



Po 



X 

2xp2Ur2 



dx{xUr2) H d4,U^2 



= 0, 



7Po 
Po 



D1P2 + 2xp2Ur2) = 0, 



X 



(28) 
(29) 

(30) 
(31) 
(32) 
(33) 



Note that Poisson's equation is no longer required, since we 
have separated out the gravitational potential perturbation 
in this approximation, which is equivalent to Cowling's ap- 
proximation ( |Cowling|1941 1. 

We can rewrite these equations in a more natural nota- 
tion, removing the asymptotic scalings to find 



contracting Eq. [42] with u 



1| ,2 

2'"' +2C^ 



+ V' 



1, ,2 & 



= 0. (46) 



is the energy density per unit 



■ po 



+ q) u is the energy flux 



Thus E 

volume, and F_b 
density. 

If the fluid is at rest, with 6 = 0, then the stratification 

surfaces are circles (spheres in 3D), and measures the 

strength of the stable stratification. If we disturb the fiuid 

from rest, then a positive (negative) buoyancy is associated 

with an inward (outward) radial displacement of particles, 

resulting in an outward (inward) acceleration of the fiuid 

due to buoyancy to restore the system to equilibrium. From 

the energy equation Eq. |46| we can see that the state 6 = is 

the state of minimum gravitational potential energy, since 

,2 

the available potential energy density poi^=^ is minimised 
for this state. This makes sense, since this corresponds to 
having a background state with no wave-like disturbance. 



DUr 



r 

UrU^ 



= —drq + rb, 
1 



pr 



-dr{rur) + -dd,Ud, = 0, 
r r 



Db + Crur 



0, 



D ^ dt + Urdr + —d^, 
r 

where 

= ^P2, 

Po 



1 

Po' 



(34) 
(35) 

(36) 
(37) 
(38) 

(39) 
(40) 



are a buoyancy variable and a modified pressure variable, 
and 



= 4<E>2 



P2_ 

IPO 



p2 

po 



(41) 



is related to the buoyancy frequency N hy N = Cr, with 
> Q for a stably stratified centre of a star. 
We now write them in the vector-invariant form 



Du = -Vg + rb, 
Db + C'^r-u = 0, 
V-u = 0, 
D = 9t + u ■ V. 



(42) 
(43) 
(44) 
(45) 



These equations are similar to the standard Boussinesq 
system for a slab of fluid in Cartesian geometry (see e.g. 
Biihler|[2009 Ch. 6) with a uniform stratification, with the 
exception that our problem is in cylindrical geometry, with 
g and A*" proportional to r. Note also that the buoyancy 
variable defined here is related, but not identical, to that 
of the standard Boussinesq approximation used in atmo- 



spheric sciences and oceanography (see e.g. Biihler 2009 1 



The buoyancy variable is proportional to the density and 
entropy perturbation. 

An energy equation for our system can be derived by 



5 LINEAR THEORY OF IGWS 

APPROACHING THE STELLAR CENTRE 

5.1 Linear solution steady in a frame rotating 
with the pattern speed of forcing 

If the radiation zone is forced from above then this will excite 
waves, which will propagate to the centre of the star. If the 
incoming wave has frequency ui and azimuthal wavenumber 
m, then it seems reasonable to assume that the response is 
steady in a frame rotating with the angular pattern speed 
of the forcing fip = a;/m, in the absence of instabilities. The 
dependence on (j> and t is then only through the combination 



(47) 



We can choose Q.p/C as a unit of length, and fip ^ as a unit 
of time (note that these units are only used in this section 
and in Appendix [A|), to allow us to write the equations in 
the dimensionless form 



DUr 



= —drq -f rb, 



Dus -\ = Oeq, 

r pr 

-dr{rur) + —deus, 
r r 

Db + rur = 0, 

D — Urdr + m ( — — 1 
V r 



(48) 
(49) 

(50) 
(51) 
(52) 



The energy equation (Eq. \46\ allows us to infer that the 
radial energy flux 

^2Tr r 1 , 

rUrd^ (53) 



po 



is independent of r for disturbances steady in this frame of 
reference, since the solutions are periodic with period 2n 
starting at ^ = 0. 

The radial angular momentum flux is 



(54) 
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We can obtain a solution to these equations by lineariza- 
tion as follows, assuming that the solution is proportional 
to e'^. Then we obtain (where real parts are assumed to be 
taken) 



-imur = —drq + rb, 
imq 
r 



1 im 

-drirUr) H Urn = 0, 

r r 
—imb + rur = 0. 



(55) 
(56) 

(57) 
(58) 



The incompressibility constraint allows us to express 
the velocity in terms of the streamfunction i/)(r, (j)), which is 
defined by 



U : 

SO we can write 
im 



Ur ~ Re 



— Re [—drip] 



(59) 

(60) 
(61) 



This enables us to reduce the system to Bessel's equation of 
order m, 



Lmip = drirdr^p) + r 1 



?/) = 



(62) 



with solution regular at the origin tp oc Jm{r). This repre- 
sents a wave that approaches from infinity, refiects perfectly 
from the centre and goes out to infinity. Pure ingoing and 
outgoing wave solutions are described by Jm{r) ± iY,n{r) 
respectively. 

5.2 Properties of the (non-)linear solution 

The general solution can be written in the form of a sum of 
ingoing and outgoing waves, with complex amplitudes Ain 
and Aout, as follows: 



V.n(r,0 = [Jm(r) + iy™(r)]e'«, 



(63) 
(64) 
(65) 



We can check that ipin corresponds to an ingoing wave 
by calculating its phase and group velocities. A simple cal- 
culation shows that the radial phase velocity is directed out- 
ward if we adopt the convention that u) = mQp > 0, and the 
group velocity is directed inward. This highlights one of the 
peculiarities of IGWs - that the phase and group velocities 
are oppositely directed, as discussed in § [2] For these linear 
waves, we have 



= pQnmrlm[ipdrip*] = 2mpo (j^ontl^ - l^inj^) 



(66) 



The solution for a wave that perfectly reflects from the cen- 
tre is 



^(r,0 = 2AnJ,n(r)e'^ 

and has — 0, since Ain = Aout- 

If we take the curl of Eq. |42| we obtain 



(67) 



9t(V X u) = V X (r6) - V X (u- Vu), (68) 
which has eliminated the modified pressure perturbation q. 



The z-component of this equation expressed in terms of the 
streamfunction is 



9i(-VV) = -d4.b + JiiP, -VV), 
and the buoyancy equation is 
dtb^ -C^d^i; + J{4;,b). 

The nonlinear terms take the form of Jacobians, 

^ diA,B) ^l d{A,B) 
d{x,y) r d(r,4>) 

= {drA){-d4,B)-{-d4,A){drB). 
r r 

The Jacobians of the solution derived above are 

j(^,-vV) = J(V',&) = o, 



(69) 
(70) 

(71) 
(72) 

(73) 



which expresses the surprising result that the solutions de- 
rived are exact nonlinear solutions of the system. This fol- 
lows from the fact that —V^tp = b — ip for these waves. This 
arises because although the nonlinear terms u- Vu 7^ 0, they 
are balanced by the modified pressure term in the equations 
of motion. We also have u ■ Vfe = 0. 

This is distinct from, but analogous to, the result that 
a single propagating plane IGW in a uniform stratification 
is a nonlinear solution of the standard Boussinesq system 



(Drazin 1977 Klostermeyer 19821. This is a consequence 



of the fact that k • u = for these waves, which implies 
that the advective operator u ■ V annihilates any distur- 
bance belonging to the same plane wave. A useful conse- 
quence of this is that a stability analysis can be performed 
on finite-amplitude propagating IGWs, allowing detailed un- 
derstanding of the initial stages of the breaking process for 
these waves ( Drazin|1977 Klostermeyer 1982 : Klostermeyer 
[1991 ). Such studies have shown that a single propagating 
IGW solution is always unstable whatever its amplitude, 
since it undergoes resonant triad interactions ( Drazin|1977 |. 
One important difference in our problem is that the nonlin- 
earity is spatially localized to the innermost wavelengths, 
whereas the nonlinearity is present everywhere in the plane 
IGW problem. We postpone a study of the stability of our 
nonlinear standing wave until a subsequent paper, though 
we expect waves of sufficiently large amplitude to be unsta- 
ble if they overturn the stratification. 

The amplitude required to overturn the stratification 



can be derived from our solution Eq. 67 The entropy (or 



more precisely, a quantity proportional to the entropy) is 
s = b + (l/2)r'^, in these units. Overturning the stratifica- 
tion means that the entropy profile, perturbed by a wave 
with buoyancy perturbation b, must satisfy drS < 0, which 
implies {l/r)drb < —1 is a condition for overturning. Since 
b = %p for these nonlinear waves, this can be expressed in 
terms of the streamfunction as {l/r)drip < —1. Reintroduc- 
ing dimensional variables, and substituting for u^, modifies 
this criterion to {u^/r) > fip, i.e., overturning occurs if the 
angular velocity of the wave exceeds the angular pattern 
speed. Equivalently, wave breaking occurs if 



max(it0) 



4C" 



(74) 



whose largest value occurs where the amplitude of drJ2{r) 
is largest, which is one wavelength from the centre. 

In 3D the equivalent linear solution (written down in 
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the Appendix of OL07), is not an exact nonlinear solution. 
However, the results of a weakly nonlinear analysis find that 
the reflection is close to perfect, with a reflection coefficient 
that is extremely close to unity for moderate amplitudes. 



6 NUMERICAL METHODS 



6.1 Snoopy spectral code 



We solve the system of equations (42 1- (45 1 using a Cartesian 
spectral code. Snoopy ( [Lesur fc Longaretti||2005j [Lesur fc] 
Longaretti|2007 L ft is advantageous to use a Cartesian code 
over one in the more natural (for the problem) cylindrical 
geometry, because of the absence of a coordinate singularity 
at the origin, near to which is the region of the flow that 
we are most interested in. We also avoid the timestep issues 
close to the centre that would be present in a time-explicit 
cylindrical code. These arise from the CFL condition, be- 
cause the grid spacing becomes very small near the origin. 

Since this is a Fourier spectral code, the problem must 
be periodic in space. We solve our non-periodic problem us- 
ing this code by setting up a region near the outer boundary, 
in which the fluid variables are smoothed to zero as we ap- 
proach the boundary, using a parabolic smoothing function. 
We find it is quite acceptable to do this over a region about 
1/10 of the total box size. For this value there is negligible in- 
teraction between neighbouring boxes. This approach is one 
that may be useful in many applications which would bene- 
fit from the use of spectral methods, but have non-cartesian 
geometry and/or non-periodic boundary conditions. The ob- 
vious drawback of such an approach is the slight increase in 
computational cost, since the smoothing region is additional 
to the flow in the region of interest. Interior to this we have 
a thin ring in which we implement a forcing term in the ra- 
dial momentum equation of the form fr cos{2(j} — tot) . This is 
designed to excite ICWs with m = 2, but is not designed to 
accurately describe the excitation of IGWs at the top of the 
radiation zone, since we are only interested in the dynam- 
ics of the central region. Our forcing is non-potential, which 
reflects the fact that the tidal forcing of waves is indirect 
(OL04; [Ogilvie||2005[ ). A potential force would be absorbed 
in this model by a hydrostatic adjustment of q. 

We solve the equations 



Du = -Vg + r6 + !/V^u+ I °' 



0, < r < rforce, 

force ^ ^ <^ smooth-) 

C^r • u = 0, (76) 
V ■ u = 0, (77) 
D = a + u ■ V, (78) 

where f = ~fr cos(2(f)—u)t) e^. We use a parabolic smoothing 
function ci(r) = r^mooth. ^ ^ , to instantaneously smooth 
UrjUtj, and 6 to zero as we approach the outer boundary. We 
do this by multiplying the variables in the region rgrnooth ^ 
r < rtox by d{r) during every timestep. 

Our choice of units for length and time are arbitrary, but 
we choose the following. We study a region —1.5 < a; < 1.5, 
— 1.5 < y < 1.5 and set rforce = 0.85 rtox, Tsmooth — 
0.9 Vhox, with Tbox = 1.5. We choose a typical IGW radial 
wavelength of \r = 0.1, so that we are resolving ~ 12 wave- 
lengths within the box. The radial wavelength of these waves 



is not strictly constant, but its variation for large r is small, 
and can be reasonably approximated by a constant value 
in that region. We then choose a forcing frequency tj = 1. 
These choices are arbitrary, and are made to ensure that we 
are resolving a sufficient number of wavelengths within the 
box. 

Explicit viscosity is added to Eq. |42| since the code 
has no intrinsic dissipation. This is necessary for stability - 
to ensure that we have no unphysical growth of energy at 
small scales. The value of the viscosity is chosen such that 
it dissipates disturbances on the grid scale, and a value of 
V — 2 X 10~® is chosen for all simulations. Viscous terms 
are implemented in a time-implicit manner. We do not in- 
clude thermal diffusion in the buoyancy equation since this 
was found to be unnecessary for stability. We solve the rele- 
vant Poisson equation for the modified pressure during each 
timestep. 

We normalise the velocity components with respect to 
a typical radial phase velocity of the wave u)\r/2TT and thus 
set 



Uj\r 



(79) 



in which is equivalent to the wave steepness s, and is 
a measure of the nonlinearity in the wave. This allows us 
to write the condition for overturning the stratification in 
Eq.[74las 



1 



(80) 



Since we have chosen to specify the radial wavelength 
Xr and frequency cj of the waves that we wish to study, we 
have already constrained the stratification 



C 



TTLO 



(81) 



There is now only one further parameter, fr (except for 
viscous damping and smoothing terms), to fully specify the 
problem. We set 



fr — fr 



LjXr 
~2^ 



(82) 



and we vary the normalised amplitude fr to model the ef- 
fects of different tidal forcing amplitudes. From preliminary 
investigation, we find that it is appropriate to choose values 
between 10~^ and 10'^, since these result in central ampli- 
tudes that range from u,- ^ 1 to ftr = 0(1) (higher central 
amplitudes are not observed, as is described in the results, 
owing to wave breaking above a critical amplitude). This 
represents a vast range of amplitudes of tidal forcing, from 
cases in which the secondary body is a low-mass planet to 
a solar-mass binary companion in a close orbit. 

Our background is a hydrostatic equilibrium with no 
wave, with 6 = in the initial state. We use a resolution 
of 512 X 512 for most simulations, though several higher 
resolution runs have been performed using 1024 x 1024 and 
1536 X 1536. We confirm that the results are not dependent 
on the numerical method (and that our system Eqs. 42|45 
correctly describes the relevant physics) by reproducing the 
basic results using ZEUS-2D (iStone & Normanll 19921). We 



describe our implementation of the problem in this code in 
Appendix [B] ZEUS reproduces the same basic results as 
the Snoopy code, which indicates that the effects of nonzero 
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compressibility are unimportant. In light of this, we only 
discuss the Snoopy results below. 



7 NUMERICAL RESULTS 



We use the set-up described in § |6.1| for a set of simulations 
with a variety of forcing amplitudes fr. The typical radial 
group velocity and wave crossing time are, respectively, 



tc = 



fbox 



(83) 
(84) 



For the initial conditions described in the previous section, 
tc ~ 90. We define for the purposes of the following, a "wave" 
to be a non-axisymmetric oscillatory flow represented by a 
single azimuthal wavenumber m 7^ 0, whereas a "mean flow" 
is an axisymmetric azimuthal flow with m — 0. 

We perform several different quantitative analyses of 
the results. We separate the amplitudes of the waves into 
an ingoing wave (IW) and an outgoing wave (OW), and 
calculate a reflection coefficient, using the method described 
in Appendix|X] The reflection coefficient TZ is defined as the 
ratio of the absolute amplitudes of the outgoing (Aout) and 
ingoing waves {Ai„) for a given radial ring. 



(85) 



and measures the amplitude decay for a wave travelling from 
r to the centre, and back to r. We can relate it to the phase 
change on refiection (A<j!>) by 



In 71 - In 



A in 



(86) 



For perfect standing waves, Ai„ = Aout, and 7?. = 1. If 
the ingoing wave is entirely absorbed at the centre, then 
7?. = 0. Thus, 7?, is a measure of how much the wave has 
been attenuated on reflection from the centre. 

We also Fourier analyse the solution, to study the tem- 
poral evolution of different azimuthal wavenumbers in the 
flow. This is done by selecting a ring of cells in the grid at a 
particular radius, which is chosen to be at r = 0.1, since this 
is probably close enough to the centre to detect the effects 
of nonlinear wave couplings, if they occur. Since this is a 
cartesian grid, we do this by selecting all cells within a par- 
ticular radial ring to within a tolerance width comparable 
with the size of a grid cell. We then compute the Discrete 
Fourier Transform of the velocity components, and from this 
calculate the power spectral density, 



p — 



1 

iV 



Ur,k exp 

fc=0 



'imk 



(87) 



and similarly for u^, where A'^ is the number of grid points 
in the ring - which depends on r and the resolution, though 
10^ < < 10^ for all resolutions at r = 0.1. From this 
we can determine which components of the solution grow or 
decay as a result of viscous damping, instabilities or nonlin- 
ear wave-wave interactions. Note that this is only a rough 
approximation to the azimuthal power spectral density be- 
cause the points are irregularly spaced around the ring, yet 




-0.1, 



Figure 2. Radial velocity along the x-axis for a simulation with 
forcing amplitude insufficient to cause breaking. The amplitude of 
this wave is largest in the centre. Also plotted is the corresponding 
analytic standing wave Eq |67| converted into a radial velocity 
using Eq |60[ showing that our simulations accurately describe the 
waves for the case in which the waves reflect coherently from the 
centre. 



we have assigned an even weighting to each point. Neverthe- 
less, this is justified in practice because this method works 
well when tested on low-amplitude solutions that are well 
described by the standing wave in § [5] for which we know 
that the solution is an m = 2 wave for both Ur and w,^. Note 
also that Pm ~ P-m since u is real, so we cannot distinguish 
between waves with wavenumbers m and — m without also 
observing the time-dependence of the fiow. 

The main result that will be discussed in more detail 
below is that we find that there exists a critical wave am- 
plitude beyond which wave breaking occurs near the centre. 
Below this amplitude, the waves reflect coherently from the 
centre, and a steady state is reached in the reference frame 
rotating with fip, consisting of an m = 2 standing wave so- 
lution. This is the outcome inferred from linear theory (see 

and we will discuss these cases first, followed by those 
in which nonlinear effects start to become important. We 
refer to the former as "low-amplitude" cases, and the latter 
as "high-amplitude" cases. 



LOW-AMPLITUDE FORCING: 
REFLECTION 



COHERENT 



When the simulations are started, transients are excited by 
the forcing at many different frequencies (and radial wave- 
lengths), centred around a; = 1 in frequency space. As more 
inward propagating waves are excited by the forcing, an in- 
going wave train propagates toward the centre. At this stage 
in the time evolution, the solution is composed of many dif- 
ferent frequencies, so our decomposition of the solution into 
a single IW and OW does not work well. As more transients 
escape the region and are damped, the primary response of 
the fiuid is in the form of waves with frequency w = 1 and 
azimuthal wavenumber m — 2. 

As the waves approach the centre and reflect, an OW 
is produced. As this process continues, the amplitude of the 
OW matches that of the IW near the centre. For these low- 
amplitude cases, the phase change on refiection is negligible. 
This means that we have coherent reflection from the centre, 
which allows standing waves to be produced. These waves 
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Figure 3. 2D plot of the radial (top) and azimutiial (bottom) 
velocity in the equatorial plane of the star for small-amplitude 
waves. This is at a time t = 36tc, once standing waves have 
formed, in a simulation with max(si0) < 0.3, in which we have 
coherent reflection from the centre. In the outer part of the grid, 
the solution is smoothed to zero to satisfy periodic boundary con- 
ditions. 



are stationary in a frame rotating with Qp, as inferred from 
linear theory. We confirm that our simulations produce the 
correct standing wave solution, by plotting an example of a 
comparison between the simulation and the wave solution 
in Fig. [2] We plot the velocity components in two dimen- 
sions for an example simulation in which standing waves 
have formed for a small-amplitude case with ma,x{u,/,) ~ 0.3 
in Fig.jS] 

After a few wave crossing times, the reflection coefficient 
increases to values approaching unity throughout the grid, 
though its value decreases with radius, shown in Fig. |4] for 
a low-amplitude case, with max({i^) ~ 0.3. In this figure 
we also plot the results of our IW/OW decomposition in a 
small-amplitude simulation with a resolution 1536 x 1536. 
Our reconstructed solutions match the data well except very 
close to the centre, thus showing that our decomposition 
works well for these cases. 







Figure 4. Radial velocity along the line y-axis (top) and az- 
imuthal velocity along the line y = x (middle top) in a small- 
amplitude simulation, after standing waves have been set up at 
t = 36tc, with max(n0) ~ 0.3 (solid lines). Also plotted is the 
reconstructed solution using Ai„ and Aout obtained using the 
method described in Appendix [A| (dashed lines). These are well 
matched everywhere except near the centre showing that our de- 
composition works well for these cases. In the bottom right we 
plot \Ai„\ (solid curve) and |Aouf| (dashed curve) versus radius, 
together with a linear fit to each curve. The waves damp as they 
propagate due to viscosity. In the bottom panel we plot the re- 
flection coefficient Ti versus radius. 



Waves and tidal dissipation in a solar-type star 11 



The decay in TZ with radius is a result of the nonzero 
viscosity, which results in a decay of wave amplitude with 
time (and therefore distance from where they are excited). 
The OW has been damped for longer, which results in the 
amplitude of the OW being smaller than that of the IW. 
A simple estimate of the amplitude decay due to viscosity 
with propagation from radius r and then reflected back to r 
again gives 

— cx exp -2 / dr Usi exp — r , (88) 

since k ~ except near the centre, and Cg^^ — cijAr/27r 
throughout the box. This roughly matches the amplitude 
decay between Ain and Aout at r = 1.2, implying that the 
decay in amplitude is indeed due to viscous damping of the 
waves. In addition, the regular oscillations in the amplitudes 
result from the fact that our exact solution in the inviscid 
case is not an exact solution in the presence of viscosity. This 
was verified by running a low-amplitude simulation with v — 
0, in which case the oscillations disappear. 

We find that the wave reflects coherently from the cen- 
tre when < 0.5. Long-term simulations {t ~ several hun- 
dred tc) do not show the development of any instabilities 
that act on waves with < 0.5, though there is a slow 
growth of m = components of in the solution, as can 
be seen in Fig. [5] In this figure, we plot P,„ for the first 
few even wavenumbers in the fiow from an example low- 
amplitude simulation. Negligible growth in odd m-values is 
observed, which is consistent with the symmetry of the basic 
wave and the quadratic nonlinearities of the Boussinesq-type 
system. The growth in m = is a result of viscosity, which 
acts to damp the waves and transfer angular momentum 
from m = 2 to the mean flow. This can distinguished from 
a process resulting from nonlinear interactions, because it 
is found to depend on i/. However, most importantly, no 
instability is observed for waves with < 0.5. 
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Figure 5. Temporal evolution of the power spectral density Pm 
at r = 0.1, in the lowest four azimuthal wavenumbers m in the 
solution for (top) and (bottom), in a low-amplitude simu- 
lation with ~ 0.3 near the centre, for grid resolution 512 X 512. 
The solution is in the form of m = 2 waves, and reaches a steady 
state in the frame rotating with Qp. Growth of m = is nonzero 
as a result of viscous damping of the waves. No wave breaking 
occurs because < 0.5 in the solution. The m = 6 components 
are most likely due to errors in the Fourier analysis. 



9 HIGH- AMPLITUDE FORCING: WAVE 
BREAKING AND CRITICAL LAYER 
FORMATION 

If we increase the value of fr, then the above picture changes 
considerably when a critical wave amplitude is exceeded. 
Once U0 > u^^crit ~ 0.5, wave breaking occurs near the 
centre within several wave periods (a few 2tt/lu), and the 
outcome of the simulations is very different from the small- 
amplitude case. This occurs when the wave overturns the 
stratification - see Eqs.[74]and[80] In Fig.[6]we plot the 2D 
velocity components after wave breaking has occurred in a 
simulation with fr = 15. 

For highly nonlinear forcing, for example with fr > 20, 
the waves break as they reach the centre with sufficient 
amplitude before there has been any significant refiection. 
For fr ~ 10, the amplitude of the IW alone is insufficient 
to cause breaking, and we must wait for refiection at the 
centre to produce an OW of comparable amplitude before 
> u^^crit- Once this critical value of is exceeded, the 
waves break. This is an irreversible deformation of the oth- 
erwise wavy material contours ( Mclntyre||2000 1 . 

Once breaking occurs, we observe consequent mean fiow 
acceleration (i.e. growth of m = components of u^), as the 



angular momentum of the waves is deposited locally where 
the wave breaks. This acts to spin up (if fip > 0) the cen- 
tral regions, which at this stage contain a sufficiently small 
fraction of the angular momentum of the star that their spin 
can be readily affected by these waves. Once this process has 
begun, the central regions spin up to ~ f2p, and a critical 
layer is formed, at which the Doppler-shifted frequency of 
the waves goes to zero. At this location, the azimuthal phase 
velocity of the waves would equal that of the local rotation 
of the fluid, if they were ever to reach it intact. In reality, 
as subsequent IWs approach the critical layer, nonlineari- 
ties dominate, and the waves undergo breaking before they 
reach it - though see discussion in § |9.2[ We plot the angular 
velocity of the fluid normalised to fip in the bottom panel 
of Fig.[7l 

As IWs approach the critical layer, their radial wave- 



length decreases, and they slow down, i.e. Cg 



as 



cj — > 0. This causes a buildup of wave energy just above the 
critical layer, in which nonlinearities become important. In 
this thin region, the quadratic nonlinearities produce higher 
wavenumber disturbances with even m values - see § |9.1[ 
These are produced by the self-nonlinearity of the primary 
IW (m = 2) as it approaches the critical layer - these self- 
nonlinearities vanish in the absence of a mean flow. These 
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Figure 6. The top panel shows the radial velocity in a high- 
amplitude simulation with /r = 15 at i = 7tc, after wave breaking 
has occurred, from a simulation with a resolution 1536 X 1536. The 
bottom panel shows at the same time in the simulation. 



daughter waves are damped faster than the primary be- 
cause they have lower frequencies and therefore shorter ra- 
dial wavelengths, which is a result of the theorem proved in 
Hasselmann ( 1967]). Thus the IW is irreversibly deformed, 
and transfers its angular momentum to either the mean flow 
or to daughter waves that are then more easily dissipated by 
viscosity. A large fraction of the angular momentum of these 
waves must be given to the mean flow when these waves dis- 
sipate. This process acts to spin up the fluid just above the 
critical layer to ~ fJp. As subsequent IWs are absorbed by 
the critical layer, the spatial extent of the mean flow expands 
outwards, i.e. the star is spun up from the inside out. We 
envisage that this process will continue until the mean flow 
encompasses the bulk of the radiation zone or the planet 
plunges into the star - though see § |10[ Long-term simula- 
tions, lasting for several hundred tc, show that this appears 
to be the case. 

This picture is analogous to |Goldreich fc Nicholson| 



(19891, who propose that early- type stars in close binaries 
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would spin down (if > n, or spin up if n < n) from the 
outside in, once a critical layer has formed near the sur- 
face as a result of radiative damping of the waves, where 



Figure 7. The top panel shows the radial velocity along y = x for 
the same case as Fig. |6] In the region r G [0.3, 1.2] the solution 
is well described by linear waves with |Aout| <C l^inli but not 
near the centre. The bottom panel shows the angular velocity of 
the fluid normalised to the angular pattern speed of the forcing 
Qp along the x-axis, at the same time in the simulation. This 
shows that the central regions after wave breaking are spun up to 
slightly exceed Qp. The critical layer occurs where u^/{rQp) = 1. 



this effect is strong. In our problem, an instability of the 
primary wave, which occurs once the wave overturns the 
stratification, causes wave breaking. This results in angular 
momentum deposition and spin up (for the case in which 
Q < n, spin down if Q > n) of the central regions, which 
causes the formation of a critical layer near the centre of 
a solar-type star. The rate of expansion of the spatial ex- 
tent of this region depends on the forcing amplitude; for 
larger amplitudes it expands faster. The critical layer moves 
outwards when the dissipation of subsequent IWs deposits 
sufficient angular momentum to spin up the fluid to ~ Qp, 
and the spatial extent of the mean flow expands. In this 
picture, there is a front of synchronization which gradually 
moves outwards (though see § |10[) . 



9.1 Growth of azimuthal wavenumbers in the flow 

We now discuss the results of a spectral analysis of the simu- 
lation data so that we can study the growth of the mean flow 
(m — 0) and daughter waves produced by breaking (other 
|m| 7^ 2 wavenumbers). We plot Pm for the first few even 
wavenumbers for a set of examples in Figs. ID and [9] Negligi- 
ble growth in odd m-values is observed, which is consistent 
with the symmetry of the basic wave and the quadratic non- 
linearities (though it is possible in principle for the wave to 
be unstable to odd-m perturbations). 

At the beginning of the simulations m = 2 dominates 
until the primary wave breaks and transfers angular momen- 
tum to the mean flow. When subsequent waves approach the 
critical layer the primary IW transfers angular momentum 
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Figure 8. Temporal evolution of the power spectral density Pm 
at r = 0.1, in the lowest four azimuthal wavenumbers m in the 
solution for (top) and ti^ (bottom), in a simulation with fr = 
15 for resolution 1536 X 1536. The solution is primarily in the 
form of m = 2 waves, until growth of even-m disturbances occurs 
as the wave breaks. Note that m = grows strongest for u^, i.e. 
angular momentum is transferred from the primary waves to the 
mean flow. Once r = 0.1 is located inside the rotating region, Ur 
is primarily composed of \m\ = 2. 

to higher m-value disturbances - this can be seen from Fig. [8] 
prior to t ~ 100, after which the ring r = 0.1 is enveloped 
by the mean flow. After this, m = dominates u^. On the 
other hand, Ur is then primarily in the form of |m| = 2 
disturbances, which from examination of simulation output, 
counter-rotate with the forcing, and have angular pattern 
speed — fip. The excitation of these waves could explain the 
counter-intuitive effect of the most central regions spinning 
slightly faster than Qp, since they carry negative angular mo- 
mentum. These waves appear to reflect from both the m = 2 
critical layer (though note that these waves do not see this 
as a critical layer) and the centre. As these waves approach 
the m — 2 critical layer, since they are counter-propagating 
waves, their frequency is Doppler-shifted upwards towards 
A'^, and they undergo total internal reflection. These waves 
appear to reflect back and forth from the critical layer and 
the centre. 

We also note the appearance of oscillations in the energy 
in m = 6 in the solution. This is most likely due to errors 
in the Fourier analysis, because we are not sampling the 
solution with evenly spaced points. The amplitude of these 
oscillations is much smaller than that of the m = 2 or m = 
waves, so should not affect any conclusions drawn from these 
results. 
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Figure 9. Temporal evolution of the power spectral density Pm 
at r = 0.1, in the lowest four azimuthal wavenumbers m in the 
solution for (top) and (bottom), in a simulation with fr = 
2.5 for resolution 512 X 512. Until t = 2200, viscous damping acts 
on the waves, transferring angular momentum to the mean flow. 
Once > u^ i-rit in the solution, wave breaking occurs. This 
occurs at t ~ 2200, resulting in a jump in the growth of m = 0, 
and a drop of energy in |m| =2. 

We experimented with the forcing amplitude to study 
cases in which wave amplitudes were just sufficient to cause 
breaking after running the simulation for ~ 25tc. The results 
of our spectral analysis of the results of such a simulation 
are plotted in Fig. [9] We clearly see evidence for viscous 
damping in producing growth of energy in m = 0. This can 
be distinguished from the sudden growth which results after 
the onset of the instability that leads to wave breaking. The 
growth of an instability on the primary wave occurs once the 
wave overturns the stratification, when U0 > U0,crit near the 
centre, at t = 2200. After onset, the wave breaks and a crit- 
ical layer is formed, leading to the general picture described 
above. For this simulation, viscous damping of the wave is 
responsible for spinning up the flow, and the formation of a 
critical layer. 

9.2 Discussion of wave reflection from the critical 
layer and implications 

Once the critical layer has formed, we find that a large frac- 
tion of the IW angular momentum is absorbed near the cen- 
tre. First, we confirmed this naively by watching animations 
of the time dependence of the velocity components. For both 
components, the wave pattern moves outwards, which corre- 
sponds to inward propagating IG Ws (see § [5| , so at least a 
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significant fraction of tfie solution is in the form of IWs. This 
is quantified by performing our IW/OW decomposition. The 
results of this for a typical simulation are plotted in Fig. |10[ 
where fr = 15. The time has been chosen after the critical 
layer has formed, and the mean flow has been accelerated 
near the centre. The IW/OW wave decomposition does not 
work well near the centre, as we might expect, since here 
the disturbance is primarily the mean flow (m = 0), though 
there are also other components. Several wavelengths from 
the critical layer, in the region 0.5 < r < 1.2, the recon- 
structed wave solution matches the simulation output quite 
well. The matching is much noisier than in Fig. |4] since the 
solution contains contributions from m 7^ 2 wavenumbera 
and uj ^ 1 frequencies in addition to m = 2, tj = 1 waves. 

Fig. [To] shows that the amplitude of the IW decays as 
it propagates towards the centre. There is significant ab- 
sorption of wave angular momentum near the critical layer, 
as the IW propagates through the mean shear. This re- 
sults in |Ao„t| <^ \A-in\, though the reflection is nonzero so 
l^outj 7^ 0. This leads to 7?. 1 in the region where the 
decomposition works well, which implies that most of the 
angular momentum in the IWs is absorbed near the centre 
- also note that the energy flux ratio oc TZ^ ^ 1. In addi- 
tion, the phase of the OW is perturbed with respect to the 
IW, which inhibits the formation of standing waves. TZ oscil- 
lates with radius, mainly because the solution is composed 
of some 7^ 1 and m 7^ 2 components, which are not filtered 
by our IW/OW decomposition. 

When a single propagating wave approaches a critical 
layer, the outcome has previously been found to depend on 
the ratio of the strength of nonlinear wave- wave couplings to 
linear viscous and radiative damping. Nonlinear wave-wave 



couplings occur over a timescale ( Booker & Bretherton 1967 1 



tNL^O k,'\drUhr^Ur 



(89) 



and linear viscous and radiative damping occur over a 
timescale 



tL=0 k7-^ 



'\drUh\ 



(90) 



The ratio of these terms define the parameter (Maslowe 
1986l|Koop|198l| 

1/3 

A = — ~ ( ^^^^ 1 . (91) 



1^ 



u\drUh\ 



where kh is the horizontal wavenumber, drUh is the typical 
shear in the mean flow, Ur is a typical radial velocity in 
the wave, and is a diffusivity - which for the centre of a 
star is likely to be primarily radiative diffusion rather than 
viscosity, so f will be primarily the thermal conductivity k. 

Nonlinearity acts to promote energy transfer away from 
the critical layer through the generation of daughter waves, 
and linear damping acts to suppress this resonant wave pro- 
duction. If A ^ 1, then the time required for nonlinear ef- 
fects to become important is long compared with that for 
linear damping to become important. In this limit, there is 
negligible wave reflection, and nearly all incoming wave en- 
ergy is absorbed by the critical layer, as predicted by IHazel] 
(19671. In the opposite limit, when A ^ 1, nonlinear ef- 







Figure 10. Radial velocity along the line y = x (top) and 
azimuthal velocity along the x-axis (middle top) in a large- 
amplitude simulation with wave breaking, at t = 7tc, with 
fr = 15. These are plotted (solid lines) together with the recon- 
structed solution using Ai„ and Aout obtained from the method 
described in Appendix [A] (dashed lines) - these are well matched 
for r £ [0.5,1.2]. Bottom right panel shows Ai„ (top line) and 
Aout (bottom line) vs radius. There is significant absorption of 
IWs near the critical layer at t ~ 0.1, resulting in |^o^if| 
The bottom panel shows the reflection coeflflcient 72. versus radius. 
72 1 in the region where the decomposition works well. 



fects become manifest prior to the time when they become 
suppressed by linear damping. In this case, the flow can 
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be extremely complicated, and nonlinearity in the critical 
layer region can lead to wave refle ction, with am pltiudes 



19711. In this 



O (exp (^-n{Ri - 1/4)^/^^^ or less (Breeding 
limit, nonlinear wave-wave couplings lead to the generation 
of many smaller-scale daughter waves, some of which prop- 
agate away from the critical layer, carrying a fraction of the 



wave energy (Fritts 19791. Experiments of internal waves 
approaching a critical layer have been performed by |Koop| 
( lisij ) and |Koop fc McGee| |l986t , in which A = 0(1). 
For this value, the effects of both nonlinearity and linear 
damping become manifest at approximately the same time. 
Their laboratory experiments show that wave reflection, in 
the form of daughter waves produced by nonlinear couplings 
that propagate away from the critical layer, is suppressed by 
viscosity for this value of A. 

Relating these results to our simulations, we find typical 
values of the parameter A ~ 10~^ near the critical layer, 
due to the explicit viscosity (u — 2 x 10~®) added in the 
code, since |9rC^h| ~ 0(10^) and Ur = Ur ~ 0.1. In this 
limit nonlinearities are likely to become important before 
viscous diffusion. Since this is the limit in which we would 
be expected to find wave reflection, if any occurs at all, and 
we find little reflection of waves from the critical layer, it is 
likely that most of the IWs are absorbed near to the centre, 
and not reflected. In any case, the reflected waves will not 
have the same frequency and horizontal wavenumber as the 
primary wave. Instead, reflected waves will be in the form of 
disturbances with smaller frequencies, and therefore shorter 
radial wavelengths, as well as higher m- values, as a result of 
wave- wave coupling ( Hasselmann|| 1967 1 . Such disturbances 
will be more easily dissipated by radiative diffusion, since 
the rate of energy dissipation ~ m^/oj*. 

We conclude that if wave breaking and critical layer for- 
mation occurs, it is probably reasonable to assume that the 
IWs are entirely absorbed in the radiation zone, primarily 
near to the critical layer. This is inferred from our simula- 
tions, as can be seen in the example in Fig. |10[ in which 
l^outl <S |j4i„|. The result of this is that if wave break- 
ing and critical layer formation occurs, it is not possible for 
global standing modes to develop in the radiation zone. It 
would then be appropriate to calculate the tidal dissipation 
rate using the method of GD98. 



10 DISCUSSION 

10.1 Orbital evolution of the planetary companion 

The main motivation for our work is to study the tidal Q' 
expected for solar-type stars, and in particular to connect 
this with the survival of close-in extrasolar planets. If the 
amplitude of the tide is sufficient to excite IGWs that ex- 
ceed the critical amplitude for wave breaking at the centre, 
then our results find it highly likely that the wave reflection 
will be imperfect, and a large fraction of the incoming wave 
energy will be absorbed. It is then reasonable to estimate 
the amount of tidal dissipation by assuming that the incom- 
ing waves are completely absorbed. This was first done by 
GD98 for a nonrotating solar model, and we here reproduce 
the relevant quantities to calculate Q' using their approach. 
The energy fiux F in ingoing IGWs excited at the convec- 
tive/radiative interface (at radius r — rint) is (see GD98 Eq. 



13) 
F 



(92) 



where oj = 2n is the frequency of IGWs excited by a planet 
in a circular, nonsynchronous orbit. 
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is a constant (for quadrupolar tides 1 = 2), and 
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where Q depends on the stellar properties at the interface 
between radiative and convective regions, and 
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dyn 
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and (Jc is a constant whose value depends primarily on the 
thickness of the convection zone, and is equal to -1.2 for the 
current Sun. The values of these quantities depend on the 
stellar model. The amplitude of the largest tide for a circular 
orbit is (e.g. OL04) 



5 {rrii, + nip) 



(97) 



if the tidal potential has the form 
Re[*r2p2™(cos6l)e'('"'^-'^*'], where corresponds to 

the Legendre polynomial of degree 2 normalised over 6 (re- 
lated to the spherical harmonic Y2,,n = (\l ^J2^)P^ e^""*). 
Converting the energy flux F to an angular momentum flux, 
and assuming all wave angular momentum is deposited in 
the star, we calculate a torque T, using the torque formula 



given in Peale (19991, giving 
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where m = 2. Then 
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(100) 



where P is the orbital period, which is twice the tidal pe- 
riod for the diurnal component of the tide. Note that the 
exact value depends on the stellar model adopted - in par- 
ticular the value of G, which is determined from the stellar 
properties at the interface between radiative and convective 
regions, and the thickness of the convection zone. The given 
value applies to a stellar model of the current Sun (see below 

for a discussion of the variation in this parameter for other 

^1 

stars). The largest uncertainty is likely to be in , 
since this may depend on the modelling of convective over- 
shoot. However, uncertainties in Q for the current Sun are 
much less than an order of magnitude, so this is a good es- 
timate of Ql that results from the mechanism described in 
this paper, for our closest star. 



16 A.J. Barker & G.I. Ogilvie 



We now estimate the orbital and stellar properties re- 
quired to excite waves that are sufBciently nonlinear near the 
centre for breaking to occur, using the Appendix of OLO'ij^ 
and the above energy flux. The nonlinearity parameter in 
OL07 is defined so that for ^ > 1, the wave overturns the 
stratification during part of its cycle. Equating the energy 
fiux in ingoing waves near the centre to Eq. [92] gives 



A = 
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where « 2 x 10'''^kg m^s^/^ and Cq 8 x 10""m-^s-\ If 
we assume that the primary instability acting on 3D waves 
occurs for the minimum amplitude at which the 3D wave 
solution overturns the stratification {A > 1), which we have 
found is indeed the case in 2D, then we have the following 
criterion for wave breaking: 



C_ 



Mj 



Mp 



1 day 



>3.3. (103) 



A Jupiter-mass planet in a one-day orbit around the 
current Sun would not raise tides of sufficient amplitude 
near the centre for breaking to occur, and would likely sur- 
vive, because it does not satisfy this criterion. Note, however, 
that Jupiter, with its period of P = 4332d, does satisfy this 
criterion. The waves excited by Jupiter in the Sun would 
be of very low amplitude, but they are also of very low 
frequency. This means that their wavelength is extremely 
short, so the energy of these waves would be concentrated 
into an extremely small volume near the centre of the star, 
if they were to reach it. However, radiative diffusion is cer- 
tain to damp these waves before they reach the centre, since 
they are of such short wavelength. This process would, in 
any case, contribute negligibly to the orbital evolution of 
Jupiter. 

We performed a study of the variation in the parameters 
Q and C using an extensive set of stellar models, with masses 
0.5 S5 m*/M0 ^ 1.1, and ages that represent the range of 
main-sequence ages expected for these stars. These were pro- 
vided by J0rgen Christensen-Dalsgaard, and were computed 
using ASTEC ( [Christensen-Dalsgaar dl|2008| ) . This involved 
integrating GD98 Eq.(3) throughout the convection zone in 
each of these models, where A'^ ~ 0, using a linear shoot- 
ing method, to determine CTc. The results of this study are 
that ~ 5© to within a factor of 5 for all solar-type stars, 
throughout the main-sequence age of each star. This is true 
even taking into account the evolution of the position of 
the interface between convection and radiation zones, and 
the resulting change in the density of the star at the inter- 
face. The main uncertainty in each model is probably | j 

at the interface, since the slope of A^ there is uncertain. 
However, this is only raised to the -1/3 power, so this effect 
probably does not contribute more than a factor of 2 un- 



^ Note that the ApJ version of this paper has a misprint in equa- 
tion (A3). The first term in the square brackets should have the 
coefficient \/6 replaced by 1/ ^6, since the given function does 
not go to zero as x — > 0. 



certainty in Q. Taken together with changes in stellar mass 
and radius, we find that our estimate of in Eq. 100 is 
quite robust, if critical layer formation occurs at the cen- 
tre. This is approximately true for all stars within the mass 
range 0.5 ^ mt/i\f0 ^ 1.1, throughout their main-sequence 
lifetime. 

The variation in A is primarily dependent on C, since 
A oc The strong dependence on C means that as 

a star evolves the tide could become nonlinear at a critical 
age, since C increases with evolution on the main-sequence 
(see Fig. [ijfor the evolution of C with the age of the Sun). 
For a given age, this parameter is found to be larger in more 
massive stars, due to their greater central condensation. In 
addition, stars with lower metallicity also have a greater 
central condensation for a given age, and so have larger C 
values over stars with higher metallicity. Over the range of 
stars considered in this study, C is found to take values 
between 0.1 — IOC©. This leads to a large variation in A 
values, for fixed orbital parameters. As a result, this param- 
eter is critical in determining whether wave breaking occurs 
at the centre. We have stated that a short-period Jupiter- 
mass planet does not satisfy Eq. |103| around the current 
Sun. However, such a planet around a similar age l.OM© 
star with a metallicity Z = 0.01 will cause wave breaking at 
the centre, since C is larger by a factor of 3. Thus, there is 
a strong dependence of the breaking criterion on the stellar 
model, primarily through the parameter C. 

Tidal dissipation of the quadrupolar tide raised in the 
star leads to evolution of the semi-major axis at the rate 
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The inspiral time for a planet into the current Sun is then 
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1 day 
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Note that the strong frequency depen- 
dence of Q'^ means that this mechanism could be very im- 
portant for short-period systems. This predicts that a planet 
spiralling into its star will undergo rapid acceleration as it 
migrates inwards. This is as a consequence not only of the 
reduction in semi-major axis, but also the decrease in 
as the tidal frequency increases with the inspiral. It must 
be noted that simple timescale estimates do not accurately 
refiect the evolution if the orbit is eccentric, inclined, or if 
the stellar spin is not much slower than the orbit ( [Barker fc| 
Qgilvie|2009 1, but this estimate shows that this mechanism 
can be very efficient in contributing to the tidal evolution of 
hot Jupiters on the tightest orbits. Indeed, Ta < 5 Gyr for a 
Jupiter-mass planet in an orbit of less than about three days 
around the current Sun, if it were to cause wave breaking 
near the centre. 

We can crudely estimate the maximum orbital period of 
a planet that can be pulled into the star by this process by 
equating the moments of inertia of the radiation zone (which 
extends to radius Rrz), to that of the orbit no? = r^girii,, 
giving 



P ~ 1.8 days 



where ^ 



Rf 



0.7 Rr. 



Mj 



(106) 



is the reduced mass and r; = 0.076 IS 
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the dimensionless radius of gyration for a polytrope of in- 
dex 3. For a planet with such an orbital period, the whole 
of the radiation zone must be spun up to cause the planet 
to completely spiral into the star. Once the entire radia- 
tion zone has spun up to the orbital frequency, this process 
becomes ineffective and the corresponding tidal torque will 
vanish. However, the tidal torque is nonzero due to dissipa- 
tion of the equilibrium tide by turbulent convection, and so 
even if this process stops, it does not guarantee the survival 
of the planet. In addition, magnetohydrodynamic coupling 
between the convection and radiation zones could act to par- 
tially counteract the spin-up of the interior, if the convection 
zone is spinning slower than the orbit. Magnetic braking of 
the star through the interaction of its magnetic field with a 
stellar wind acts to spin down the convection zone, which 
would also gradually spin down the radiation zone through 
these couplings. If the coupling between the convection and 
radiation zones of the star is efficient, then our mechanism 
would again become effective. In the case that the coupling 
timescale balances the timescale for spin up of the radiation 
zone due to IGW absorption at a critical layer, then the 
planet would migrate into the star on the magnetic braking 
timescale. A detailed study of these effects is not currently 
possible, since there are many uncertainties, but this is wor- 
thy of future consideration. 

A Jupiter-mass planet on a ~ 1 day orbit will spin up 
a substantial fraction of the radiation zone on infall. If the 
ratio of the orbital moment of inertia to the spin moment 
of inertia of the radiation zone > 1, then the above process 
alone will be unable to cause the planet to spiral into the 
star. 



10.2 Critical layer formation induced by radiative 
diffusion 

Near the centre of a star, radiative diffusion is the dominant 
linear dissipation mechanism. If the waves are sufficiently at- 
tenuated by radiative diffusion on their reflection from the 
centre, then over a sufficiently long time, they may be able 
to spin up the innermost regions of the star to the angu- 
lar pattern speed of the tide Sip, hence producing a critical 
layer. We have already observed this process occuring in our 
simulations in Fig. |9] albeit with viscosity and not radiative 
diffusion acting on the waves. Here we provide an order-of- 
magnitude estimate of the timescale for this process. 

If we assume that the wave is attenuated as it propa- 
gates from a radius R, to the centre, and back to R again, 
by a factor e~" , i.e., that the energy flux is Fe~°', then the 
torque on this region of the star is given by the angular 
momentum transferred to the mean flow, and is 



UJ 



(107) 



Here we assume that F is that given in Eq. |92[ which is 
probably reasonable when the response is non-resonant. The 
wave attenuation a is given by 



2 r "^^^dr. 



(108) 



where Jjrad is the thermal conductivity, k is the wavenumber. 
We have k ~ kr except within the last wavelength from the 
centre, and the radial group velocity Cg^r ~ The thermal 



conductivity can be calculated from the properties of the 
appropriate stellar model from 
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where a is the Stefan-Boltzmann constant, k is the opacity, 
and Cp is the specific heat at constant pressure. We can take 
rjrad to be constant over the inner ~ 3% of a star, to a first 
approximation, with a value « 16.7m^s~^ in the case of the 
Sun. In addition, k,. is roughly constant with radius within 
this region, though diverges as r — > 0. We can reasonably 
estimate 
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where R is the size of the region spun up by this process. 

if < AT, and N = Cr. We 



This is because kr 

also note that the tidal frequency oj — 2 (^). This means 
that the attenuation by radiative diffusion within the inner- 
most few percent of the star, is small for waves excited by 
planets on one-day orbits. However, the wavelength of the 
waves becomes shorter for longer period orbits, so their at- 
tenuation by radiative diffusion is more efficient. Note that 
a > 1 only when P > 8 days, even when radiative diffusion 
over the entire radiation zone is considered (in fact GD98, 
who made a more accurate calculation by including the ra- 
dial dependence of the wavenumber and diffusivity, found 
that a > 1 for P > 11.6 days). This means that the waves 
excited by planets on short-period orbits with P < 3 days, 
whose survival could be threatened by the process of wave 
breaking, will not be significantly attenuated in traversing 
the radiation zone, according to linear theory. It is therefore 
appropriate to ask whether they will break on reaching the 
stellar centre. 

To a first approximation, the central ~ 3% of a star 
can be modelled as a uniform density sphere, with the cen- 
tral density pc. Its moment of inertia is 7 = ^ivR^pc. Using 



Eq. |107| we have the following differential equation describ- 
ing the spin evolution of the central regions: 

7^ = -P(l-e-"). (Ill) 
at UJ 

Assuming that the system evolves slowly, which is proba- 
bly true until breaking occurs, we can take F, lj, and a, 
to be constant in time as the central regions are spun up. 
This allows a straightforward solution, giving the resulting 
timescale to spin up the region of the central wavelength to 
Qp, in the case of Jupiter orbiting the Sun with a one-day 
period (in which R ~ O.O17?0), of 



4 ttR^PcQ.1 

tsu = ^- 
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3 Myr. 



(112) 



Note that this timescale strongly depends on the size of the 
region that is spun up. When a <C 1, which is valid for 
P< 10 days, tsu oc P^^^. 

In this estimate we have used Model S of the current 
Sun, so this value only applies to our star. However, we have 
repeated this calculation using the same set of stellar models 
discussed in the previous section, with masses in the range 
0.5 m*/MQ ^ 1.1, and find that tsu ^ 1 Gyr for each of 
these models, for Jupiter orbiting the star with a one-day 
period. 
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This striking estimate indicates that all gas giants on 
short-period orbits around G or K stars could eventually 
cause the formation of a critical layer near the centre of the 
star, given sufficient time < 0(1) Gyr. Once this has formed, 
we have found in our simulations in i|7j that it is reasonable 
to assume that the ingoing wave angular momentum flux is 
entirely absorbed near the centre. Hence, our estimate of Q'^ 
in the previous section could apply to all slowly rotating G 
and K stars. The effects of rotation will complicate matters 
if the tidal frequency is less than twice the spin frequency, 
as studied in OL07. 

Unlike the mechanism of nonlinear wave breaking, 
which is the main subject of this paper, the formation of 
a critical layer by radiative diffusion requires the progres- 
sive spin-up of the region of the central wavelength by a 
very gradual deposition of angular momentum. This pro- 
cess could be interrupted by other mechanisms of angular 
momentum transport that resist the development of differen- 
tial rotation, such as hydrodynamic instabilities or magnetic 
stresses. A consideration of such effects will be required in 
order to determine whether this process operates in reality. 



10.3 Long-term evolution of the radiation zone 

As the planet migrates inwards due to the IGWs it excites 
being absorbed at a critical layer, cj increases, since n > 0. 
If we assume that the tidal frequency increases from ui at 
a time ti to 102 > oji at a slightly later time t2, then since 
ct;2 > 1^1, the waves excited at t2 no longer see the critical 
layer for lji waves. If the change in frequency is sufficiently 
small, we would still expect significant attenuation by the 
shear as the waves approach the critical layer for uji waves. 
This would transfer angular momentum from the waves to 
the mean flow, and may spin up the region near the origi- 
nal critical layer to the pattern speed for 102 waves, hence 
producing a critical layer for these waves. If this process oc- 
curs, then a weak radial differential rotation profile could be 



set up in the radiation zone, with 



dn(r) 



> 0. Similarly, if 



the planet migrates outwards, then a profile with 



< 



could be set up. If the flow does not have time to adjust to 
the change in frequency of the forcing, and cannot spin up 
sufficiently to produce a critical layer for UJ2 waves, then the 
dissipation rate may be reduced. However, it seems plausible 
that the change in the orbit will be gradual enough so that 
IGWs reaching the centre will be significantly attenuated. 
This is because their radial wavelengths get Doppler-shifted 
downwards by the shear, making them more susceptible to 
radiative damping. 

We have so far restricted our investigation to m = 2 
waves. If the orbit of the planet is eccentric or inclined with 
respect to the stellar equator, then IGWs with other m- 
values could be excited. If these break, then they could each 
have their own critical layer. It would be interesting to study 
the effects of tidal forcing with several different frequencies 
and m-values to see how this would affect the reflection of 
the different waves. In addition, 'Rogers & Glatzmaier ( 2006 1 



found that IGWs excited by turbulent convection at the top 
of the radiation zone, which had 1 < m < 15 and typical fre- 
quencies u) ~ 20/^Hz, could reach the centre with sufficient 
amplitude to undergo breaking and spin up the central re- 
gions. If this commonly occurs in solar-type stars, then we 



might expect the core to be differentially rotating, even in 
the absence of tidal forcing. For most m-values the result- 
ing spin frequency of the central regions is probably slower 
than that of the relevant pattern speed for tides raised by a 
close-in planet, so we expect that this would minimally at- 
tenuate any tidally excited low-m waves. Nevertheless, the 
interactions of multiple waves near the centre warrants fur- 
ther study. 



11 CONCLUSIONS 

We have presented a study of the fate of internal gravity 
waves approaching the centre of a solar-type star, primar- 
ily using two-dimensional numerical simulations. A train of 
internal gravity waves excited at the interface between the 
convection and radiation zones propagates towards the cen- 
tre. These waves break if they reach the centre with steep- 
ness sufficient to overturn the stratification, which in 2D 
corresponds to > u^^crit ~ 0.5. Once this occurs, nonlin- 
ear effects cause the subsequent formation of a critical layer, 
as the waves transfer their angular momentum to the mean 
flow, bringing the central regions of the star into corotation 
with the tidal forcing. This acts as an absorbing barrier for 
subsequent ingoing waves, which continue to be absorbed 
near the critical layer, resulting in an expansion of the spa- 
tial extent of the mean flow. The general picture of this 
process is that the star spins up (or down) from the inside 
out, until either the planet has plunged into the star, or the 
radiation zone of the star has spun up (or down) to match 
that of the evolving orbit. 

In 2D, wave breaking occurs near the centre if > 0.5. 
The amplitude of gravity waves excited by tidal forcing have 
this value near the centre of the current Sun if 



Mj 



1 day 



> 3.3. 



(113) 



However, this value depends strongly on the stellar model, 
in particular the value of ^ at the centre, which varies 
between different stars, and with main-sequence age. It is 
likely that a short-period Jupiter-mass planet will not ex- 
cite waves with sufficient amplitude to cause breaking near 
the centre of the current Sun, since it does not satisfy this 
criterion. There is also a dependence on the stellar properties 
at the interface between convection and radiation zones (see 



Eq. 103 1. However, the uncertainties in these parameters are 



expected to be markedly less than an order of magnitude, 
given a stellar mass and approximate age (and metallicity) . 

By decomposing the numerical solutions into an ingoing 
and outgoing wave, we find that if wave breaking occurs most 
of the angular momentum of the ingoing wave is absorbed 
near the centre, and is not reflected. This has very important 
implications for tidal dissipation in solar-type stars. We find 
that the assumption of GD98 and OL07, that tidally ex- 
cited internal gravity waves approaching the centre a star 
are not coherently refiected, and do not produce standing 
modes, is appropriate. Neglecting the effects of rotation, we 
use the calculations of GD98 (with the minor corrections 
of OL07) to estimate the modified tidal quality factor that 
results from dissipation of these waves near the centre. Its 
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value is 

Q[ ^ 1.5 X 10^ 



1 day 



(114) 



for the current Sun. This mechanism produces efficient tidal 
dissipation over a continuous range of tidal frequencies, once 
wave breaking occurs. If wave breaking does not occur, the 
wave is perfectly reflected from the centre, and global stand- 
ing modes can be set up in the radiation zone. In this case, 
efficient dissipation only occurs when the tidal frequency be- 
comes resonant with a global standing mode, but this con- 
tributes negligibly to the dissipation, since the system moves 
rapidly through these resonances ( |Terquem et al.||1998[ ). 

From studying an extensive set of stellar models, this 
estimate of is found to vary by no more than a factor of 
5 between all main-sequence stars, with masses in the range 
0.5 — I.IM0, at any stage in their main-sequence lifetime. 
This means that our estimate is quite robust, and is likely 
to apply to all G and K stars within the mass range 0.5 — 
I.IM0, at any stage in their main-sequence lifetime (as long 
as they do not possess a convective core). 

The strong frequency dependence of Q'^ implies a rapid 
(~ 1 Myr) and accelerating inward migration of planets on 
the tightest orbits around solar-type stars, if these planets 
excite waves with sufficient amplitudes to cause breaking. 
These planets can spiral completely into their stars if their 
orbital moment of inertia is smaller than the spin moment 
of inertia of the radiation zone. If the planet has a larger 
moment of inertia than that of the radiation zone, then the 
planet will spin up the entire radiation zone, and this pro- 
cess will become ineffective in the absence of any compet- 
ing effects. Consequently, we predict that fewer hot Jupiters 
which satisfy the breaking criterion, with orbital periods of 
less than 2-3 days, will be found around solar-type stars. 
Planets with masses much larger than Jupiter could trans- 
fer sufficient angular momentum to the radiation zone before 
they are engulfed by the star. Further evolution would then 
depend on the strength of the coupling between the radia- 
tion and convection zones, and on magnetic braking of the 
star through the interaction of its magnetic field with a stel- 
lar wind. The outcome that results from the interaction of 
these effects is uncertain. 

Our work is an extension of nonlinear mechanism for 
tidal dissipation proposed by GD98, and contributes several 
new results. While their model was applied to the circular- 
isation of close binary stars, we (following OL07) focus on 
the inward tidal migration of short-period extrasolar plan- 
ets. We have performed numerical simulations, which clearly 
determine the criterion for wave breaking in two dimensions. 
Furthermore, we have found that the resulting spin-up of the 
central regions of the star leads to an ongoing tidal torque 
as incoming waves are absorbed in a critical layer that mi- 
grates outwards. Our work allows a calculation of the tidal 
quality factors in an extensive range of stellar models when 
this process occurs. It also raises issues for future investiga- 
tion, particularly concerning processes of angular momen- 
tum transport in the radiation zone and their interaction 
with the wave breaking process. 

Of the very close-in hot Jupiters currently observed, 
those around W ASP- 12 ( |Hebb et al.|2008t, WASP-18 pleF 
|lier et al.|[2009t and OGLE-TR-56-b ( |Sasselo^[2003t , the 
host stars all have convective cores. This means that the 



mechanism outlined in this paper will not work in such stars, 
regardless of the mass of the planet. This is because the con- 
vective cores would reflect internal gravity waves and pre- 
vent them from reaching the centre, where they can break. 
For such stars, the tidal dissipation is found to be weak 
( Barker fc Ogilvie|2009 1, potentially explaining the survival 
of those planets. 



A recently discussed hot Jupiter is CoRoT-2b (Gillon 
et al.|2009 1, for which Spitzer IR observations have been re- 
ported. This system has a nip ~ 3.5Mj planet in a P = 1.74 
d orbit around a star with a radiative core, with mass 
m* ~ 0.96Mq. The age of this star has been inferred to 
be young, due to the high Li I abundance, its rapid rota- 
tion, and the strong Ca II H & K lines, though there are 
significant uncertainties. However, best estimates give the 
main-sequence age < 300 Myr. At this young age, we ex- 
pect the nonlinearity of the gravity wav es ap proaching the 



centre of the star to be weak (see Eq. 101 1, since at 



the centre is likely to be small compared with that of the 
current Sun (see Fig. [T]). This means that the process that 
we have discussed is unable to operate, and based on these 
considerations, we would expect Q'^ to be large, providing a 
potential explanation for the survival of this planet. 

In addition, a recently accepted paper reported the dis- 
covery of the latest candidate for the shortest period tran- 



siting planet, around the G-type star WASP-19 (Hebb et al, 
[2010| ). This planet has mass rup = 1.15Mj, in an orbit of 
P = 0.78 d, around a star of mass ~ 0.95Mq, and 
hence will contain a radiative core. The stellar age is poorly 
constrained, but we find that the gravity waves excited by 
this planet will not have sufficient amplitude to cause wave 
breaking at the centre of the star, for all reasonably aged 
stellar models of a similar mass star. This means that the 
planet will not be subject to accelerating tidal decay through 
this process, perhaps explaining its survival. Constraints on 
the tidal Qi for this system would then give us information 
on alternative mechanisms of tidal dissipation, such as the 
dissipation of the equilibrium tide by turbulent convection. 

We have identified a further process by which efficient 
tidal dissipation of internal gravity waves could occur. If the 
waves are of too low amplitude to initiate breaking, the weak 
deposition of angular momentum through radiative damping 
of the waves, can spin up the region of the central wavelength 
over a timescale of millions of years, until a critical layer is 
formed and the mechanism discussed previously can con- 
tinue. This process could provide efficient tidal dissipation 
in solar-type stars perturbed by less massive companions. 
However, it may be prevented by hydrodynamic or magne- 
torotational instabilities, or by other effects that resist the 
development of differential rotation near the centre of the 
star. 

This work has pointed out the importance of nonlinear 
effects near the centre of a solar-type star in contributing to 
tidal dissipation. Linear theories of tides (e.g. OL04; OL07; 
[Savonije fc Witte|[2002l |Witte fc Savonije||2002| must in- 
clude a correct parametrisation of the effects of internal wave 
breaking near the centre (like that in OL07 and GD98), if 
they are to correctly determine the efficiency of tidal dissi- 
pation in solar-type stars. So far, this work has been only 
two-dimensional, has included only a single component of 
the tidal potential, neglected any rotation profile in the ra- 
diation zone, as well as any possible influence of magnetic 
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fields. In future work, we plan to study the stability of the 
nonlinear gravity wave solution in § |5] to understand the 
initial breaking process in more detail. In addition, three- 
dimensional effects could be studied in an extension of our 
current simulations. This appears to be a promising mecha- 
nism of tidal dissipation, which warrants further study. 
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APPENDIX A: REFLECTION COEFFICIENT 
CALCULATION 

Since we have an exact analytic solution to the Boussinesq- 
type problem (see !|5|, we can can deconstruct the numerical 
solution into a sum of ingoing and outgoing waves. Doing 
this enables us to quantify the amount of angular momentum 
absorbed as these waves approach and/or reflect from the 
centre. The approach we use is now described. 

A solution can be decomposed into the sum of an in- 
going and outgoing wave as (in the dimensionless units of 

t/'in('-,C) = 

V'oi,t(r, = 
V.(r,C) = 

where ^ = m4 
are 



[Jm{r) + iYm{r)y^ , (Al) 
[Jm{r) - iYm{r)y^ , (A2) 
Re [Air^ipin {r, + Aoutipout {t, 0] , (A3) 
-Lot. The velocity corresponding components 



u,.(r,C) 



= Re 



— At„il>i„[r,^) H AoutiPout{r,0 

r r 



(A4) 



Re[vrir,0 + Wrir,0], (A5) 
Re [-Aindripin [v, (_) - Aoutdri^out {v, ^)] ( A6) 
ReK(r,C) + w</,(r-,e)], (A7) 



for general m £ Z. From now on we choose m — 2, which 
corresponds with that of our forcing. These are steady, but 
non-axisymmetric, in the frame rotating with Jlp = oj/m. 
We assume that the amplitudes are locally independent of 
r, so that we can ignore their radial derivatives. Since the 
amplitudes Ai„ and Aout are in general complex, we require 
4 pieces of information from the simulations to be able to 
calculate them. We choose to calculate these from Ur and 
given from the simulations at two different azimuthal posi- 
tions, around the same radial ring. We repeat this process 
for 8 pairs of points around the ring, and take the mean of 
these values to calculate Ain and Aout- This involves solving 



where x = (Ain, Aout) (which has 4 com- 
ponents for complex amplitudes) and b — 

{ur(r, ^i) , U4,{r, ^i) , Ur(r, ^2) , U4,(r, ^2)) , and the matrix 



A 



Ro[i,r(r, £1)] 
Reb^(,-,5l)l 



Im[i,r(r, 5i)] 



Ra[™^(r,ii)] 
Re[™^(r-,«l)] 



Im[™^(r,{i)] 



Re[«r(^, «2)1 «2)] Re[™^(r,{2)l Im[™* (V,' {2 ) j ) (^9) 



Ro[«^(r, 52) 



»["^('-, 52)1 Ro[™^(r, 52)1 



52)1 



for any two points with (p — and = ^2. If we choose 
^2 — ^1 = for n £ Z, this matrix is non-singular for 

all radii, and is related to the Wronskian. This is calculated 
for all radial rings in the grid. We verifled this method on an 



analytic standing wave solution (Eq 67 1 using Mathematica, 
and wrote a Matlab code to read in Snoopy/ZEUS data and 
compute Ain, Aout and the reflection coefficient TZ, which is 
defined in Eq. |85[ For perfect standing waves, Ain — Aout, 
and TZ — 1. If the ingoing wave is entirely absorbed at the 
centre, then TZ = 0. 

The main reason for this approach is that it enables 
us to compute all four unknowns (the real and imaginary 
parts of the complex amplitudes) for each r, which requires 
Ur and for two azimuths. Integrating azimuthally over a 
ring to eliminate m 7^ 2 components, would leave only two 
unknowns, Ur and - which is not sufficient to determine 
the amplitudes. The disadvantage of our approach is that 
m 7^ 2 components also contribute to the amplitudes. This 
problem can be ignored if we trust the computed values of 
TZ, only where the solution is well described by the linear 
solution i.e. far from the wave breaking and forcing regions. 

When the solution has different frequency and 
wavenumber components than a; = 1, m = 2, the solution 
cannot be simply decomposed into an ingoing and an outgo- 
ing m — 2 wave solution with an amplitude that is roughly 
independent of radius. In which case, the amplitudes of the 
waves that are calculated from this method if the solution is 
not predominantly m = 2 can oscillate wildly with radius, 
because the solution is not a simple m = 2 wave solution. 
If the frequency (and hence radial wavelength) of the wave 
is different from that of the chosen ingoing/outgoing wave 
solution, then we will also not be able to match the solution 
exactly. One way of gauging how well our decomposition into 
a single ingoing/outgoing wave is to plot the solution recon- 
structed from the calculated Ai„ , Aout against the numerical 
solution from Snoopy/ZEUS output. If these are very differ- 
ent then we cannot make definitive conclusions about the 
solution from these values. 



APPENDIX B: ZEUS COMPARISON 

We confirm that the results are not dependent on the numer- 
ical method by reproducing t he basic results using a s tripped 
down version of ZEUS-2EP] |Stone fc Norman|l992[ ). ZEUS 
solves the equations of ideal compressible hydrodynamics, 
using a simple Eulerian method based on finite-differences, 
implemented using a covariant formalism, enabling the use 
of non-cartesian orthogonal coordinate systems. For our 
problem we solve the problem using cylindrical (r, tf)) coor- 
dinates, which are the most natural to use for our problem. 
However, the coordinate singularity at the origin requires 



Ax = b 



(A8) 



which has kindly been made freely available by J.Stone at 



http://www.astro.princeton.edu/~jstone/zeus.html 
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that we cut out a small region at the centre, on which we 
impose reflecting boundary conditions. What may seem an 
advantage of this coordinate system, that is the higher reso- 
lution at the centre, which is automatically present when we 
use a uniform grid in r and (f), requires very short tirncsteps 
when the resolution is increased, as a result of the CFL sta- 
bility constraint. This becomes prohibitive as we increase the 
resolution of the grid to above 100 x 150 in r, (f> respectively, 
so only preliminary low resolution runs were performed using 
this code. This is also because this code solves the compress- 
ible equations, and therefore resolves sound waves. For our 
problem we require a ratio of sound speed to radial group 
velocity of gravity waves x ~ Cs/cg^r ~ 6 x 10'^ {x^^ is a 
measure of the importance of effects of compressibility), in 
order to reproduce an equivalent set-up to that used in the 
Snoopy code above, so most of the computational time is 
spent resolving sound waves. In the paper we only analyse 
the Snoopy results, since they arc at a much higher reso- 
lution, but here describe our problem set-up in ZEUS, for 
completeness. We use a circularly symmetric parabolic den- 
sity stratification, p(r) = po — p-ir^ and calculate the pres- 
sure (p) profile from hydrostatic equilibrium. We solve the 
equations 

Dp = -pV • u 

^ I* 0, Tinner ^ T* <C T forces 

Du=--'Vp+g+< f, Vforce ^ r < Vdamp, 

z,(£).-^v.„, 

where e = (7 — l)p is the specific internal energy of the 
gas, f = —fr cos{24> — cot) Br, and d = — d(r)u. We choose 
7 = 5/3, as appropriate for a monatomic ideal gas. Radial 
gravity has been implemented as a source term in the radial 
momentum equation as g = —girer- Both the inner and 
outer boundaries have reflecting boundary conditions, and 
we also implement a linear frictional damping in a region 
adjacent to the boundary to prevent the reflection of (most 
of) the outgoing wave energy from the outer boundary. A 

parabolic smoothing function d(r) — I *" ^d°"ip \ used 

\ ^box ^ damp / 

in the damping terms. We choose rbox = 1-0, rdamp = 0.9, 
r force = 0.85 and ri„ner = 0.01. In the code we specify 
pa, p2,uj,\r &i fr\ the Other relevant parameters are calcu- 
lated from these. Choosing po = 1-0, P2 = 0.1, a; = 1.0, Ar = 
0.1 and a suitable value for is sufficient to fully specify 
the problem. 

A minimum value of x = 6285 is found from these ini- 
tial conditions. Such a high value is required for the wave- 
length of the gravity waves to be — 0.1, which allows 
~ 8 wavelengths to be resolved within the grid. This value 
is not much smaller than that appropriate at the centre of a 
solar-type (x ~ 10'' — 10®). We set up the initial conditions 
in such a way to minimise this value given the above input 
parameters. 

Calculations were performed in an incrtial frame, 
though the results were interpreted in a frame rotating with 
the angular pattern speed of the tide, f2p = uj/2. In this ro- 
tating frame, the linear wave solution is steady, which allows 
the instability to be easily recognised as departures from a 
steady state. 

With this resolution, there are some numerical errors 



near the inner boundary. This results from the fact that we 
only remove a small region near the centre, which is compa- 
rable with the size of a grid cell. In addition, the code has no 
explicit viscosity or thermal conduction, so we have less con- 
trol over the scales of dissipation, than in Snoopy. Neverthe- 
less, ZEUS reproduces the same basic results as the Snoopy 
code, which indicates both that the effects of nonzero com- 
pressibility are not important in this problem, and that our 
basic results are not dependent on the numerical method. 
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